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We study the partially asymmetric exclusion process with open boundaries. We 
generalise the matrix approach previously used to solve the special case of total asym- 
Q\ . metry and derive exact expressions for the partition sum and currents valid for all 

values of the asymmetry parameter q. Due to the relationship between the matrix al- 
gebra and the g-deformed quantum harmonic oscillator algebra we find that g-Hermite 
polynomials, along with their orthogonality properties and generating functions, are 
of great utility. We employ two distinct sets of g-Hermite polynomials, one for q < 1 
and the other for q > 1. It turns out that these correspond to two distinct regimes: the 
previously studied case of forward bias (q < 1) and the regime of reverse bias (q > 1) 
where the boundaries support a current opposite in direction to the bulk bias. For 
the forward bias case we confirm the previously proposed phase diagram whereas the 
case of reverse bias produces a new phase in which the current decreases exponentially 
with system size. 
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1 Introduction 

The asymmetric simple exclusion process (ASEP) [|IJ is a much studied model both from 
a mathematical and physical viewpoint. The model comprises particles hopping in a pre- 
ferred direction on a lattice with hard core exclusion imposed. In the mathematical liter- 
ature the interest lies in its being a simple realisation of interacting Markov processes @ 
and much progress has been made in proving existence theorems, invariant measures and 
hydrodynamic limits |l|]. Early applications concerned biophysical problems such as single- 
filing constraint in transport across membranes PJ and the kinetics of biopolymerisation M. 
More recently the ASEP has achieved the status of a fundamental non-equilibrium model 
due to its intimate relation to growth phenomena and the KPZ equation ||, the problem 
of directed polymers in a random media and its use as a microscopic model for driven 
diffusive systems and shock formation |§ . Last but not least many traffic flow models are 
based on variants of the ASEP |1 . Adding to its appeal is the fact that many exact results 
have been obtained, particularly in one dimension, allowing an analytical understanding of 
the non-equilibrium phenomena exhibited [JT0| . 



In recent years the open boundary problem, where particles attempt to enter at the 
left of a one-dimensional lattice of N sites with rate a, hop to the right under a hard core 
exclusion constraint and exit at the right with rate /3, has been of considerable interest. 



It was first pointed out by Krug that boundary-induced phase transitions can occur [ ! I 



and it has further been shown that spontaneous symmetry breaking takes place when two 



oppositely moving particle species are introduced [12 



The door was opened to the analytical study of the open boundary model and its non- 
equilibrium steady state in [jnj. There, exact recursion relations on the steady state weights 
were obtained and density profiles and currents were worked out exactly for the case a = 
(3 = 1. Also a mean field phase diagram in the a- (3 plane was derived. Then the recursion 
relations were used to calculate correlation functions for the case a = (3 = 1 |TJ]]. In [15 



a different method, to be referred to as the matrix approach, was proposed: it was shown 
that the steady state weights can written as a product of matrices which are in general of 
infinite dimension and non-commuting. The matrices obey algebraic rules which replace the 
recursion relations found in Jl3| . This approach gives the full solution of the model including 



the phase diagram and density profiles and allows, in principle, calculation of all equal time 
correlation functions. Further it admits generalisation to other models. It should be noted 
that the phase diagram and density profiles were also obtained independently working from 
the recursion relations [[U|. Subsequently the matrix approach has been used to solve new 
models with several species [[T7|, [18|, [19], ^(J and various updating schemes J21], ^, |23|, to 



calculate shock profiles on an infinite system [^J] and to recover some previously known 
results 



The phase diagram obtained for the one species open boundary case |TB|, [TS], |T(J comprises 
three phases: a low density phase, a high density phase and a maximal current phase where 
generic long range correlations occur. This phase diagram appears quite robust for driven 
diffusive systems with open boundaries |26[]. Also it has been shown that different types of 



updates the same generic phase diagram pertains, except in simple cases of deterministic 



bulk dynamics where the maximal current phase is absent [27l ES 



The partially asymmetric exclusion process is a generalisation of the model where par- 
ticles are able to hop to the left as well as to the right. The case of partial asymmetry is 
of interest since it allows one to interpolate between symmetric exclusion which can have 
equilibrium steady states and the far from equilibrium asymmetric system. In the context 
of growth phenomena the two different systems are described by Edwards- Wilkinson and 
KPZ universality classes [^TJ and the crossover phenomena has been of interest [ 29j [30 1. 



A parameter q gives the ratio of hopping rates to the left and to the right; thus q = 
recovers the totally asymmetric process and q = 1 gives the fully symmetric case. Further- 
more in the range < q < 1 there exists what we shall refer to as the forward bias regime 
where particles hop preferentially to the right and the boundary conditions are such that a 
steady state current of particles to the right is supported. In contrast, when q > 1 particles 
can only enter at the left and leave at the right but hop preferentially to the left in the bulk. 
A steady state current of particles to the right can then be supported by the boundary 
conditions against the bulk bias. This gives rise to a new phase which we analyse for the 
first time here. We show that the current decreases exponentially with the length of the 
system. This phase is of interest in the context of 'backbend dynamics' |3T, 32] where, for 



example, a fluid in a permeable medium has to traverse a pore oriented against the direction 
of gravity. It also appears that such a reverse bias phase is relevant to recently reported 
one-dimensional phase separation [|3"5|, [T5 . 



In JT5) the generalisation of the matrix approach to the partially asymmetric exclusion 



process was pointed out. This was used in |34] to obtain approximate expressions for the 



current in the forward bias case for large system size. It was also pointed out that the 
matrices for the partially asymmetric case are closely related to creation and annihilation 
operators of the g-deformed harmonic oscillator. In |35| the quadratic algebra was studied 
and curves in the parameter space were deduced for which finite matrices could be used. 
Along these curves exact expressions for physical quantities, such as the current and corre- 
lation length, can be computed easily. Also, in the symmetric case q = 1 the steady state 
is straightforward to solve 0, |37|, |38[. However exact expressions for all system sizes and 
parameters have remained elusive. 

In this work we build on the relationship of the quadratic algebra with the g-deformed 
harmonic oscillator algebra to calculate exact expressions for the current for all system sizes 
and parameter values. The basic approach is to use known properties, such as orthogonal- 
ity and generating functions, of the g-deformed Hermite polynomials associated with the 
quadratic algebra. 



This method was employed independently by [39| to study the case of forward bias 
(q < 1) and a phase diagram for this regime was obtained by examining the asymptotics 
of the normalisation (partition sum). We have been able to consider for the first time 
the case q > 1 for which a less well-known set of g-Hermite polynomials is required []4U|J . 
Furthermore we obtain as a main new result of our work an exact, explicit expression for 
the normalisation valid for all q which encompasses all regimes: q = (total asymmetry), 
q < 1 (forward bias), q = 1 (symmetric) and q > 1 (reverse bias). This exact expression 



allows all physical quantities of the model to be evaluated for any system size and we use it 
here to obtain the asymptotic form of the current in the reverse bias regime. 

The paper is organised as follows: in section 2 we define the model we consider and in 
section 3 we review the matrix approach and the related quadratic algebra. In section 4 we 
discuss the g-deformed harmonic oscillator and its relation to the quadratic algebra of the 
present problem. In particular, we present relevant facts such as the generating functions 
and orthogonality relations for the g-Hermite polynomials. In section 5 we derive our main 
results, which are exact expressions for the normalisation (partition sum). In section 5T we 
derive an integral expression valid for q < 1 fl4l|) and in section |5.2| we derive an integral 
expression valid for q > 1 (f43|). We then obtain in section |5.3| a finite sum expression valid 
for all q (|57|). In section || we use the exact expressions to calculate the asymptotic behaviour 
of the current and the phase diagram. We conclude in section [7| with a discussion. 



2 Model definition 

The microscopic dynamics of the model are specified by four rates at which certain events 
can occur. For a rate A associated with a particular event, the probability that the event 
happens in an infinitesimal time interval dt is Adt. Furthermore, moves that would lead 
to two particles simultaneously occupying a single lattice site are prohibited due to the 
hard-core repulsion between them. 

The events defined in the model and the rates at which they take place are as follows. 

Event Rate 

Particle inserted onto the left boundary site a 

Particle removed from the right boundary site j3 
Particle hops by one site to the right 1 

Particle hops by one site to the left q 

Figure [l] shows two typical particle configurations on a small lattice along with the 
allowed moves and their rates. 

0^\ A A A 1 A A A ^ 



Figure 1: Typical particle configurations and allowed moves in the model. 

As only three of the rates are independent, we have set the right hopping rate to 1 with 
no loss of generality in the following analysis. 

Later, in section [1], we will consider a more general parameter space where particles can 
also enter at the right and exit at the left. 



3 The matrix product formulation and its quadratic 
algebra 

In this section we review the matrix approach to finding the steady state of the model. 
We present here the bare essentials of the method and refer the reader elsewhere for more 



detailed descriptions of the technique fTR [10], 1 1 



Consider first a configuration of particles C and its steady-state probability P(C). We 
use as an ansatz for P(C) an ordered product of matrices X\X 2 . . . Xn where Xi — D if site 
i is occupied and X{ = E if it is empty. To obtain a probability (a scalar value) from this 
matrix product, we employ two vectors (W\ and \V) in the following way: 

Pi c) = <"Wr,...*«|v) . (1) 

Z N 

The factor Z^ is included to ensure that P(C) is properly normalised. This latter quantity, 
analogous to a partition function, has the following simple matrix expression through which 
a new matrix C is defined: 

Z N = (W\{D + E) N \V) = (W\C N \V) . (2) 

Note that if D and E do not commute P{C) is a function of both the number and 
position of particles on the lattice, as expected for a non-trivial steady state. The algebraic 
properties of the matrices can be deduced from the Master equation for the process [|l^]. It 
can be shown that sufficient conditions for equation (|J) to hold are 

DE-qED = D + E (3) 

a(W\E = (W\ (4) 

(3D\V) = \V). (5) 

One can also write expressions for ensemble-averaged quantities in terms of matrix prod- 
ucts. For example the current of particles J through the bond between sites % and i + 1 is 
given by 

, _ (W\C l -\DE - qED)^- 1 - 1 ^) Z N -r 

7 ~ 7 \ ' 

where the last equality follows from relation (|3]). We see that, as expected in the steady 
state, the current is independent of the bond chosen. Also, the mean occupation number 
(density) of site i may be written as 

(W\C i - l DC N - t \V) 
Z N 

Our task now is to evaluate the matrix products in the above expressions for Zjy, J and r» 
by applying the rules 



In [[L5| the case q = was treated by using (|3|) repeatedly to 'normal-order' matrix 
products: that is, to obtain an equivalent sum of products in which all E matrices appear 
to the left of any D matrices. Then finding a scalar value from (|4]) and (|5]) is straightforward. 
For example one can develop (^) as 

Z N = (W\C N \V) = Y,aN,n,m(W\E n D m \V) = (W\V)J2a N ,n,ma~ n ^ m . (8) 

n,m n,m 

The difficulty with this approach lies in the combinatorial problem of finding the coefficients 



o.N,n,m- An alternative approach proposed in [15 is to find an explicit representation of C and 



decompose the vectors (W\ and \V) onto the eigenbasis of C to evaluate the normalisation. 

In the present work we employ mainly the latter approach to derive an expression for 
the normalisation which is valid over a restricted range of the model parameters. We will 
later compare this with the canonical form (|S|) to find an expression which does not depend 
on the chosen representation and is therefore generally valid. The representation of C that 
we use is intimately related to the g-oscillator algebra and the eigenbasis of C is constructed 
from g-Hermite polynomials. 



4 The g-deformed harmonic oscillator and its relevance 

Much progress can be made if the algebra of the previous section is written in terms of that 
of g-deformed quantum harmonic oscillator B2J . The relationship central to this algebra is 



the g-deformed commutator 

aa) — qa)a = 1 (9) 

where the operators a and a) operate on basis vectors \n) (with n — 0, 1, 2, . . . ) as follows: 

at| n ) = (t — 1 — | n + i) (io) 



l-q 

i 

a\n) = [ ~ q ) \n-l) (11) 



i-q 

a\0) = 0. (12) 

In terms of these new operators, the matrices introduced in the previous section can be 
written as 

D = -^— + ^^=a (13) 

l-q v/T^ l ; 

E = J— + —^tf (14) 

l-q y/T=q 

and one finds using (§) that @ is satisfied. 



We now have an explicit representation of the original D and E matrices in the oscillator's 
"energy" eigenbasis {|n}}. Thus one may combine equations (|TBD, QTTJ) and (|5|) to find the 
corresponding representation of the vector \V): 

W) = un "L—- (is) 

where we have set (0|V) = 1 and v is the following combination of the model parameters: 

!-9 -, 







(16) 



The denominator of (|T|) is more conveniently written in terms of a 'g-shifted factorial' 
This is defined through 



n-l 



(a;q) n = l[(l-aqi) (17) 

(a;q)o = I- (18) 

We will later encounter products of these factorials for which we shall use the shorthand 
employed by |f43fl : 

(a, b, c, . . . ; q) n = (a; q) n (b; q) n (c; q) n ... . (19) 

We can now use this compact notation to write expressions for both (n\V) and (W\n), 
the route to the latter being the same as that to find (n\V): 

(n\V) = -f= (20) 



where v is given by (|16|) and 



(W\n) = -== (21) 

Vw;<?)n 



w = l -^- - 1 . (22) 

a 



We see that the representation of D and E ( |T3| ) and (O) breaks down for certain choices 
of the model parameters. Firstly, for q = 1 (symmetric exclusion) a number of singularities 
appear. Secondly, if v > 1 and q < 1, the vector element (n\V) fl2Tf ) is unbounded from 
above as n -^ oo; similarly with (VT|n) (^1|) when w > 1 and q < 1. We consider for 
the moment only those regions of parameter space where this representation converges, and 
discuss the generalisation to the remaining areas in section |5.3|. 



We persevere with the relationship between the original quadratic algebra and the q- 
oscillator algebra for the reason we now explain. We introduced earlier a matrix C which 



appears in the expressions for the mean particle density and current. We now see that 
this matrix can be written as a linear combination of the identity 1 and the "co-ordinate" 
operator x = a + a): 

C = D + E = - 1 + - x . (23) 

1-g jT^q 

The eigenstates of the oscillator in the co-ordinate representation are known — in analogy 
with the solutions of the undeformed oscillator they are called the continuous g-Hermite 
polynomials j|3] . Clearly, the eigenvectors of C are the same as those for x and therefore 



knowledge of them permits diagonalisation of C. As this is a major step towards obtaining 
the exact solution of the model, it is worth spending a little time discussing the g-Hermite 
polynomials. 

The recursion relation for the polynomials follows after a suitable definition of the op- 
erator x on its eigenbasis {|x)}: 

x\x) = \x) . (24) 

From this and equations (|TT|) and ([K]) we find 

2x(x\n) = a/1 - q n (x\n - 1) + y/l - g n+1 (x\n + 1) . (25) 

Explicit formula? for (x\n) can be found using a generating function technique, the details 
of which differ slightly depending on whether q < 1 or q > 1. Here we present the results 
which will be most useful later; derivations are given in appendix A. 

When q < 1 we make a change of variable 

x = cos6. (26) 

The g-Hermite polynomials can now be defined as (0\n), that is, the projection of the 
oscillator energy eigenstate \n) onto the position basis (9\. The generating function 

A" 






G(M) = ^ -==<0|n) (27) 

„=o yyTiQjn 

can be expressed as an infinite product 

G ^ = v^^z (28) 

when |A| < 1. An explicit form of the g-Hermite polynomial (0\n) can be determined from 
the generating function and is presented in appendix A, equation (|A6| ) . It can also be shown 
p3| that the set of g-Hermite polynomials are orthogonal with respect to a weight function 
i/(0). That is 

r d6(n\6}v(6)(6\m}=8 n , m (29) 

o 



where 

(a e 2W e~ 2ie -a) 
v(d) = {q ' ' ^ . (30) 

Similar results emerge when q > 1 under a different change of variable 

x = isinhw (31) 

with a suitably re-defined generating function valid for all A: 

G(u, A) = V Mn) = (zg^Ae", -^Ae""; g-% . (32) 

„=0 V W ?)* 

An explicit expression for (u\n) is given in ( [A10Q . Again a weight function v{u) that or- 
thogonalises the polynomials can be found [Q. We write 



du (n\u)v{u)(u\rn) = 5 n ^ m (33) 

where now 

1 1 

V[U) ~ LT^ (g-l,-g-l e 2-,-g-le- 2 -g- 1 ) 00 ' ( } 

5 Exact expressions for the normalisation Zjv 

5.1 Integral representation for q < 1 

The relationships in the previous section allow us to obtain integral representations of matrix 
products. Here we illustrate how to apply the procedure to obtain an expression for the 
normalisation Zn when q < 1. The procedure used to obtain this particular result is also 



described by [39]; in the next subsection we extend the method to the case q > 1. 



First we take the orthogonality relation for the g-Hermite polynomials (|29| ) and use it 
to form a representation of the identity matrix: 

r 

d0|0) 1/(0) (0| = 1L. (35) 

o 

We now insert this into the expression for the normalisation (|2|): 

p'K 

Z N = d9 1/(0) (W\C N \9) (9\V) . (36) 

Jo 

By design, the matrix C is acting on its eigenvectors, so using ( p^) and (p6|) we obtain 

Z N = £mv{9)(W\6) (^l±}lj {6 \v) . (37) 



It is necessary to decompose the boundary vectors (W\ and \V) onto the {\9)} basis. By 
inserting a complete set of the basis vectors {\n)} we find 

oo oo n 

(0\V) = j>|n)<n|V) = J2^f=T m ■ (38) 

„=0 n=0 VWi I)™ 

The final sum in this equation is just the generating function of the g-Hermite polynomials 
(p7|). Thus, when \v\ < 1, we may write 

m) =G(0,v)= iveie ^ q)oo - (39) 

Similarly, when \w\ < 1 we find 

w) ^ G( ,, m) ^ __L 5 _ . (40) 

Putting all this together, we arrive at an exact integral form for the normalisation 

N p7T 

N , 



Z N =[-—\ / d6v(6)[2(l + cos 6)fG(6,w)G(6,v) (41) 

which, written out more fully, reads 

(q;qUf 1 \ N /■*,•■ , miiv (e**, e" 3 "; g) e 



^iV 



. , d9[2(l + cos 9)f -, .„ v ' J - v • (42) 



When |f| > 1 or |u>| > 1 equation (|4"T|) is not well-defined because G(9,X) does not 
converge when |A| > 1. Rather than finding a representation of the quadratic algebra that 
does not suffer from this problem, one can simply analytically continue the integral (|4~2|) to 
obtain Z^ when |u| or |iw| takes on a value greater than one. This procedure is carried out 
in section || 



5.2 Integral representation for q > 1 

We now apply the procedure of the previous section to find an integral representation of 
the normalisation for the case of q > 1 which has not previously been considered. The 
only difference is that we must use the q > 1 forms of the generating function and weight 
function of the g-Hermite polynomials, namely equations (|32H34f). We obtain a similar form 

for Z N to (0): 



(——) J duv(u)[2(l + ismhu)] N G(u,w)G(u,v) . (43) 



( i 

J N 
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However the full form is somewhat different from 
1 1 / 2 ^ 



\nq(q 1 ;q ^ \1 - q / 

[°°dv (l + i sinh v) N fa'^e", -iq-He- u , iq- l we u , -iq- l we~ u - <T% 
loo M j (-q-^-q-ie-^q-iU " ^ 



One should note the range of integration is infinite and that, in contrast to equation 
it cannot be simply replaced by a closed contour in the complex plane. It is this feature 
which makes this integral unsuited to approximation using the saddle-point method as we 
discuss section 1ST 



5.3 Explicit formula 

In this section we derive an alternative expression for Z^ which takes the form of a finite 
sum rather than an integral and is valid for all values of the model parameters. Such an 
expression is useful for two main reasons: firstly it allows us to extract the asymptotic form 
of the normalisation when q > 1; secondly, as the sum contains a finite number of terms, it 
can be evaluated exactly by numerical means if one wishes to study finite-sized systems. 

We will work from the integral for q < 1, \v\ < 1 and \w\ < 1 fl4T|). The first stage of the 
calculation is to state an important identity: 

dfl v(6)G(6, X)G(9, v)G(6, w) = - — *- — . (45) 

(vw, Xv, Xw; g)oo 

We do not prove this here but note that it is in fact a special case of the Askey- Wilson 
g-beta integral jKJ. f\ 



We now expand both sides of (j45j ) in powers of A. We already know the expansion of 
the left-hand side because G(9, X) is the generating function of the g-Hermite polynomials. 
The right-hand side may be treated by using another important identity valid when |x| < 1, 



q< 1 |3 






We find 



oo .„ n 

= y A y 

(Av, Aw; 4)00 ^ (q;q)nf^ Q 
where we have used the g-binomial coefficient which is 

(9; q)n 



(q;q)n-k(q;q)k 



v n ~ k w k (47) 



(48) 



5 In its most general form the Askey- Wilson g-beta integral has four parameters, whereas identity 
has only three. 

11 



when < k < n and zero otherwise. In the limit q —* 1 the g-binomial coefficient is equal 
to the conventional version m familiar from combinatorics. Thus the last summation in 
(f!7|) may be considered a g-deformation of the binomial expansion of (v + w) n . We give this 
function the symbol B n (v, w; q): 



B n (v,w;q) = J^ 



fc=0 



v n k w k 



(49) 



If we now compare co-efficients of powers of A on both sides of ( f4"5|) we obtain a key 
result: 



1 B n (v,w;q) 



d9u(9)(9\n)G(9,v)G(9,w) = , , 

(^;g)oo ^(q;q) n 



(50) 



This relationship is important because we may take any sufficiently well-behaved function 
f(6), re-express it it as a sum of g-Hermite polynomials and use (|50|) to evaluate the integral 



d9v{9)f{9)G{9,v)G{9,w) . 



Specifically, we can choose f{6) = [2(1 + cos^)]^ and solve equation ( f4J"l) exactly. 



(51) 



Expanding the cosine function in this way involves little more than routine algebra which 
is detailed in appendix B. The identity which emerges is 



N 



[2(1 + cos6)} N = J2 R nAi) VimU (0\n) 



n=0 



with 



I N-n I 

rnM = E (- 1 )" 

fc=0 



2N 



N -n-2k 
which may be alternatively written as 



n + k — 1 
fc-1 



+ q K 



J '/ 



n + k 
k 






fc=0 



2N 
N - n - 2k 



2N 
N - n - 2k - 2 



gl 2 J 



n + k 
k 



(52) 



(53) 



(54) 



We may now insert the expansion fl52|) into ([y]) and integrate using (p0|) : 

.AT AT 



'JV 



i / i V 



(55) 
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This exact formula, valid for q < 1, \v\ < 1 and |iu| < 1 admits extension to general q, v 
and w. We first note that the infinite product in the prefactor can be replaced with (W|V), 
a fact which follows from ((461): 



{vw;q) c 



E 

n=0 



[VW] 



(q;q)t 



(W\V) . 



(56) 



We claim that the resulting expression for the normalisation 

.AT N 



'JV 



<W|V) 



1 



1-q 



y^ j R Ntn (q)B n (v,w;q) 



(57) 



n=0 



and B n (v,w;q) by fl49|) , holds for all choices of the model 



where RN,n(q) is given by (| 
parameters. 

This is justified by observing that once v and w are written in terms of a and (3 using 



(|16|) and ( p2|) we obtain a power series in a and /? of the same form as (|8|). As discussed in 
section |3], an equation with this structure arises when one re-orders a matrix product directly 
using the relation (|3]). To perform this direct manipulation, it is not necessary to employ a 
specific representation of the quadratic algebra. Therefore, although the representation we 
used to derive ( |57| ) breaks down for q — 1, q < 1, \v\ > 1 or q < 1, \w\ > 1, we can now say 
that had we used one which converges in the region of interest, we would still have obtained 
equation ( |57| ) for the normalisation. 



As a check of this formula, let us consider the case q = 0. Then 

RN,n(0) = 

B n (v,w;0) = 



and 



2N 
N -n 



2N 
N -n- 



,n+l 



W 



ra+1 



V — W 



become 

(58) 
(59) 



where now v = 1/(3 — 1 and w = 1/a — 1. It can be verified, using the identity 

X-Y 

Em: 



X 

Y -n 



-1 



X-l-r 
Y-r 



(60) 



that ( p?D can be rewritten as 



N r 



Z N = (W\V)J2 



k=0 



2N - 2 - k 

N-k 



2N 
N- 



-2-k 
2-k 



:m 



fe+i 



:v«) 



fc+i 



l/p-l/a 



(61) 



which is equivalent to equation (39) of |15| . 



We note that for q — > 1 the singularity in the denominator of ( |57D is cancelled by the 
sum over R N ^ n (q)B n (v,w;q) and the expression is in fact well-behaved. Although we have 
checked for small system sizes that ( |5"TD agrees with the expression of [37] we have not been 
able to show this in a simple way. 
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6 The phase diagram of the model 

We now have enough information to obtain an exact phase diagram for the model and 
expressions for the particle current in the large system size limit. The behaviour differs 
greatly according to whether the particles are forward biased (q < 1) or reverse biased 
(q > 1) and so we treat the two cases separately. 



6.1 The forward bias regime 



When q < 1, the quantities of interest are most quickly obtained from the integral (42) as 
was also done in p5| . As N becomes large, we can use the saddle-point method to evaluate 
Zjv, and so we re- write ([42]) as a contour integral 

.AT 



'N 



(g;g)c 



1 



q 



dz 



(2 + z + z 



-1\N 



\z%z 



q)c 



K 



(vz, wz, vz~ 1 ,wz~ 1 ; q) c 



(62) 



where the contour K is the circle \z\ = 1 and is directed anti-clockwise. Furthermore it 
passes through the saddle-point of 2 + z + l/z along the path of steepest descent. We find 
from the saddle-point formula that 



^N 



(q;q) 



3 

oo 



^k (v,w; q) 
which holds as long as \v\ < 1 and \w\ < 1. 



3 

1 N? 
N 



1-q 



N 



(63) 



We can treat other values of v and w using the integral fl62]) after realising that (|6q ) gives 
the contribution from the contour K for any value of v and w, whereas the analytic contin- 
uation of Zn is obtained by distorting the contour such that the poles at z — v , qv, q 2 v, . . . 
and z = w,qw, q 2 w, ... stay inside it and all other poles are outside it. Figure || illustrates 
how to modify the contour as v is increased to a value 1 < v < 1/q. 




Figure 2: The integral over the contour K (shown dotted in the right-hand figure) can be 
determined from a saddle-point expansion, whereas the distorted contour K' is the correct 
one to use when 1 < v < 1/q. Note that only the four poles closest to z = 1 have been 
shown for clarity. 
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The difference between the two results can be calculated using the residue theorem. 
When \w\ < 1, 1 < v < 1/q we find by looking at the poles at z — v and z = 1/v that we 
should add 



(v 2 ;?)oo 



2 + v + v- 1 \ N 



, . , , (64) 

to (p5|). In the large N limit, this correction dominates the contribution from the contour 



K and so we write 



-2. 



«)c 



2 + w 



.-lytf 



Zw ~ 7 / — r - — i 65 

when |u>| < 1 and 1 < v < 1/q. As f is increased above 1/g and other poles need to be 
considered, (0) remains the dominant contribution to Z^. One could guess this from the 
fact that the pole at z — v is furthest from the origin. 



Due to the symmetry of (p2\) in v and w, we obtain (|65|) with v *-* w when w > 1, \v\ < 1. 
When both w and v are greater than one, the leading term in the asymptotic expansion 
comes, as before, from the pole furthest along the real axis of the complex plane. Thus we 
have found three different forms for Zjy, each of which corresponds to a phase in the model. 
These forms and their regions of validity are presented in table 0. 



Region 



v < 1, w < 1 



v > w, v > 1 



W > V , w > 1 



Normalisation Z 



N 



4 (q;q)lo (1 



7r 0,^;<?);L \ N J \ l ~q 



(v 2 ;q) 



2 + v + V 1 



N 



N 



(v w, w/v; ?)«, V 1-9 

(w" 2 ;g)oo f2 + w + w~ r '' 



(wv,v/w; q) c 



1-q 



Table 1 : The normalisation for different values of v and w when q < 1 

We can go on to find the currents in the three forward bias phases through equation (|3|). 
These expressions are presented in table 0. For completeness we should determine Z^ along 
each of the phase boundaries (v — w > 1 or v — l,w ^ 1 etc). We find that the currents 
subsequently found are equal to the limiting values of those in table |2] as the boundary under 
consideration is approached from each of the neighbouring regions. 

We can now draw a phase diagram for the system when q < 1: see figure |3]. We note 
that it has the same structure as that found for q = [15|]. The region \v\ < 1, |iu| < 1 is a 
maximal current phase; the remaining two phases correspond to the high and low density 
phases found by Jl5| . Each of these two latter phases may be subdivided into three regions 
according to the behaviour of the density correlation length in the thermodynamic limit 
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Region 



Current J 



a < i^ 2 , (3 > a 



(3< i=i, a> (3 



1 


-9 

4 




a(l- 


-g- 


a) 


1 


-q 




0(1- 


-q- 


0) 



1-q 



Table 2: The N — * oo forms of the particle current in the forward biased phases. 




Figure 3: The phase diagram of the model when q < 1. The thick solid line is a first order 



transition and the thin solid lines are second order transitions in the sense of |I5| , [16|j 



6.2 The reverse bias regime 



We turn now to the case q > 1. We will show shortly by examining the exact formula ([57]) 
that the normalisation behaves like Z^ ~ q± N for large N . This sheds further light on why 
the saddle-point method is not applicable here: by its nature, it gives expressions where the 
exponent is linear in iV rather than the desired quadratic. 

To proceed we must find approximate forms of two elementary quantities: the g-shifted 
factorial and the g-binomial coefficient. We rewrite the definition ( |TTD as 



(q;q) n = (-ir q ( n t 1 )e M ^ 



(66) 



where 



M, 



oo 
fc=l H 



- q~ kn ) 



(67) 
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We see that when n is large, M n (q) can be approximated by 

1 1 
k q k — 1 



1 1 



fc=i 
which is independent of n. This leads to an approximation of the g-binomial coefficient 



qk{ n-k) e -M( q ) ^ 



which is valid when both n and k are large. The rest of the analysis is not particularly 
illuminating, so is presented in appendix C. We ultimately find that when q > 1 



where 



Z N ~A(v,wiq)(q- 1 vw,l/vwiq- 1 ) 00 [l—j) q*" (70) 



A(v,w;q)^J^e W ^M(q) + ^^} . (71) 



We should note that the smallest system size N for which (|70D holds will be a function of 
q due to difference between the exact quantity M n (q) and the approximate form we used 
M{q). 

The expression for the current follows quickly from ([|). This reads 

which, in contrast with the currents in the forward bias regime, is a function of the number 
of lattice sites N. 



7 Discussion 

In this work we have employed properties of g-Hermite polynomials to calculate exact steady 
state properties of the partially asymmetric exclusion process. The connection between this 
model and g-Hermite polynomials lies in the fact that the matrices D and E of the matrix 
product solution can be written in terms of g-raising and q- lowering operators as in fll3|). 
Then the calculation of the normalisation @j amounts to decomposing the vectors (W\ and 
\V) onto the eigenvectors of the matrix C which are the eigenvectors of the "co-ordinate" 
operator of the q-deformed oscillator. This allowed us to obtain integral representations of 



the normalisation for both the forward bias case q < 1 ([11]) and the reverse bias case q > 1 
(fEf). Further we could use orthogonality properties of the g-Hermite polynomials to express 
these two integral expressions as a finite sum valid for all q (|57D. 
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In a very recent paper [p9|] , orthogonal polynomials were also used to study the ASEP. 



In view of the fact that our work and [B^] were carried out independently, a comparison is 



in order. In ||39| , the large system size limits of the normalisation and the current in the 
forward bias regime q < 1 via the integral representation pT| ) were obtained. Furthermore 
it was shown that one could also analyse the density correlations from this integral. In the 
present work we have found that a corresponding integral can also be found for the case of 
reverse bias q > 1. Further for all values of q (and a, (3) we have succeeded in obtaining an 
exact sum formula valid for all system sizes. One application of this general expression was 
to calculate the current in the reverse bias phase. 



For the forward bias case the phase diagram proposed by Sandow |34j] is recovered. In 
that work the more general parameter space including rates 7 (exit of particles at the left 
boundary) and 5 (entry at the right) was considered. For this case the algebra is modified 
0to 

DE-qED = D + E (73) 

(W\(aE-jD) = (W\ (74) 

((3D-5E)\V) = \V). (75) 

In principle we can generalise our method to that case. The only difference being that (W\ 
and \V) specified by @ and @ will have more complicated expressions than ( pot ) and (|2l"l) . 
In the forward bias case the generalisation would not produce any new phases. However for 
q > 1 allowing particles to exit at the left and enter at the right (both 7, 5 > 0) would allow 
a left flowing current of particles to be sustained, thus destroying the reverse bias phase. 

The reverse bias phase where the boundary conditions impose a current opposite to the 
bulk bias realises a new phase in the ASEP where the current decreases exponentially with 
system size. A typical arrangement of particles in this phase is a lattice full at the left end 
and empty at the right end. The form for the current j ~ q~ N ' 2 (|72|) suggests that the 
lattice is typically half-full i.e. the furthest particle to the right has typically to traverse a 
distance of N/2 sites against the bias to exit the lattice. To understand why the lattice is 
half full one invokes the particle hole symmetry that the current of particles exiting to the 
right must equal the current of holes exiting to the left. The symmetry implies that the 
lattice must be half- full [31 . 

Further, for large N the current tends to zero and one can compare with the much simpler 
case where the current is exactly zero, for example when the boundaries are reflecting [[TJ]. 
Then the microscopic dynamics obey detailed balance and the unnormalised probability of 
a configuration of M particles at positions x%,X2, ■ ■ ■ ,xm will be proportional to q~^^ lX \ 



Here also the normalisation grows exponentially in N [33 



In a recent preprint [[45] the density profiles for the three forward bias q < 1 phases 
were calculated in the thermodynamic limit. In particular in the maximal current phase 
it was shown that the density profile decays with distance x from the left boundary as 
1/2 + (47rx) -1 ' 2 for large x. There remain, however, a number of issues to be resolved. As 
q tends to 1 the maximal current phase occupies more and more of the phase diagram — see 
figure [| However, at q = 1 we know that the profile is exactly linear. This implies a non- 



trivial limit q — > 1 and therefore non-trivial crossover phenomena from the asymmetric to 
the symmetric case. This corresponds to the transition between KPZ and EW universality 
classes in the related growth models. Further in the reverse bias case (q > 1), as we expect 
the lattice to be roughly half-full, the density profile should be sigmoid-like. For q — > oo 
the profile will be a sharp step function whereas as q — > 1 the sigmoid profile will straighten 
out into a linear profile. 



It is known that the quadratic algebra (0-|g) for the open boundary problem can be 
used to solve a partially asymmetric periodic system with the addition of defect particles 
fl7| , |10|1 . A defect particle hops forward with rate a but is overtaken (and moved back a site) 
by normal particles with rate (3. When q < 1, the different phases in the present problem 



manifest themselves in this model [46]. However the case q > 1 is yet to be tackled; we 



believe this to be of special interest as the reverse bias phase corresponds to phase separation 



into pure domains and spontaneous breaking of translational invariance [33, O 
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Appendix A: Generating functions of the g-Hermite poly- 
nomials 

In this appendix we explain how to obtain the generating functions G(6, A) and G(u, A) 
from the recursion relations for the g-Hermite polynomials Q25|). We also present explicit 
expressions for the polynomials as they may be obtained easily from the generating functions. 

We consider first a general form of G, suitable for both q < 1 and q > 1: 

G(x,A) = ^-==(x|n). (Al) 

n=0 VW9J« 

We now obtain a functional relation for G(x, A) by multiplying both sides of equation (p5|) 
by A n /(y/(q; q) n ) and performing the required summations: 

G{x,q\) = (A 2 -2Ax + l)G(x,A) . (A2) 

By using this relation repeatedly, we can find an expression for G(x, A) in terms of G(x, 0). 
This latter quantity is fixed by normalisation, and so we set it to 1. 
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It can be seen from (|A2|) that our approach from G(x, A) to G(x, 0) depends on whether 
q < 1 or q > 1. Consider first the case q < 1. It is useful to make the change of variable 
x = cos 6 so that (|A2| ) becomes 

G(e,q\) 



G(e,\) 



(1 - Xe ie )(l - Xe- ie ) 



Iterating this we find 



G(6,\) 



1 



(A3) 



(A4) 



(Xe ie ,Xe- ie ;q) c 

where we have used G(9, 0) = 1. The infinite product l/(x;g)oo has a well-known series 
representation [E3] valid for a; < 1, g < 1 



(x;g) c 



E 

n=0 



X 



(9;?)r 



(A5) 



from which we may extract the form of (0\n). Expanding both sides of (^4) in A and 
comparing co-efficients we find 



(0\n) 



E 



3 «(n-2fc)6» 



J '/ 



V(qTq)n k=0 

where [?] is the g-deformed binomial described in section [O 



(A6) 



The case q > 1 proceeds in the same way. We must however make a different change 
of variable x = i sinh u because otherwise (|24|) would imply that we had found imaginary 
eigenvalues of a Hermitian matrix. Also we should divide A by q in as we approach G(u, 0) 
from G(u, A). We thus re- write ( |A2| ) as 

G{u, A) = (1 - iq~ l \e a ){l + iq- l \e~ u )G{u, q^X) (A7) 

and iterate before to obtain 

G(u, A) = {iq~ l Xe u , -iq^Xe^; g" 1 )^ . (A8) 

Again the infinite product on the right-hand side of this equation has a useful series expan- 
sion appropriate for q > 1 and all x: 



(q X x;q % = J^ 



x 



n=0 



(q;q)r 



(A9) 



Expansion of XH in powers of A and comparison with the generating function fl32|) yields 
(u\n): 



(u\n) 



ViqTqU 



E(-d 



fc=0 



3 (n— 2k)u 



(A10) 



J y 



It is important to realise that the two forms (|A6|) and (|A10|) we have found are not 
very different. In particular one can obtain the form ( |A10| ) by making the substitution 
9 —> 7r/2 — iu which is another way of describing the replacement cos# — > isinhu. Also 
we should note that all the functions we have found are real on their domains despite the 
presence of i when q > 1. 
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Appendix B: Expansion of the cosine function in q- 
Hermite polynomials 

Here we show how to re-write [2(1 + cos6)] N as a sum of g-Hermite polynomials (0\n). First 
of all we simplify the task by using an identity easily verified by induction: 

N f ?N \ 
[2(l + cos9)] N = J2( N _ n )c n (9) (Bl) 

n=0 ^ ' 

with 

c*(0) = \\ , m n = n (B2) 

We now need only to consider the expansion of c n (9). It is fairly easy to convince oneself 
by inspecting ( |/I6| ) that only those g-Hermite polynomials of the same parity as the cosine 
function c n will appear in the expansion. Also we do not expect any contributions from 
(9\k) with k > n. This leads to the following prescription: 

HI 

c n (6)= y %2a n , k (9\n-2k). (B3) 

fc=0 

A formula for a Hj k may be found by applying the orthogonality property of the g-Hermite 
polynomials (^9T). We obtain a familiar integral transform 

a n , k = / d6is(6)c n (6)(n-2k\6) . (B4) 

Jo 

To evaluate this integral we need the series expansion of the weight function v{6) 

1 oo 1 oo 

up) = ^yi h**® c 1 + ? s ) e " e = ?- E & ^ • ( B5 ) 



27T ^ ' 7 v 2tt 

s=— oo s=— oo 



Inserting this and the explicit formula for the g-Hermite polynomial ( |A6| ) into the above we 
find, after some manipulation, 



oo n—2k r 



®n,k — — / / b s y 

V (g; g)n-2fc s= _ tX3 r=0 



n-2£; 



2- 



— / d9 cos((k + r - s)9) . (B6) 

, 2 ^ Jo 



The integral that appears in this equation is just a representation of the Kronecker delta 
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symbol 5 Sy k+ r - Thus we can eliminate the summation over s: 



n-2k 



a nk = ( ~ 1} V 

V (?; q)n-2k r=0 



■ Ifqh) 



, k (k\ n-2k 



E 



n-2A; 
r 

n-2A; 
r 



ffc + r\ 



-l) r gV ^(1 + g 



k+r\ 



-iy q (l)+ kr (1 + /+-) 



J y 



(-l)ty 



((g fc ;g)n-2fc + g fe (g fe+1 ;g)n- 2 fc) • 



>/(?; q)n-2k 

To do the last step, we have made used the following series expansion in x: 



(z; g) n = 5^ 



n 



fc=0 L J 9 



{-l) k q®x k . 



(B7) 



(B8) 



The latest expression allows a little simplification by noting that 

{a; q)n-k 



(aq \q) n - 



2k 



and so we find we find 



a n ,k = {-l) k \J{q\q)i 



-2k q K 



{a;q)h 



n — k — 1 
fc-1 



+ q" 



J '/ 



n — k 
k 



A it 



(B9) 



(BIO) 



which is true for all n, k if we observe the usual convention that [£] =0 when k < or 
k > n. 



We may now combine equations (|B1|), ( p3|) and ( B10| ) to find 



[2(1 + cos( 



liV 



- / 2/V \ ^ 

n=0 ^ ' fc=0 



I N-n I 

AT L^~J 

EE 



n=0 fc=0 

A' 



2iV 



N - (n + 2fc) 

J^.Rjv,n(g)\/(g;g)n(%) 



On+2fe,fe(6 l |«) 



(Bll) 



n=0 



thus completing the derivation of equation (|52|). The problem of expanding of a general 
function in g-Hermite polynomials was first solved by Rogers in 1894. His approach, detailed 
by |47]], is more complicated than ours as he did not use the orthogonality properties of the 
g-Hermite polynomials. 
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Appendix C: Approximation of the normalisation for 
large N and q > 1 



We indicate here how to estimate Z^ as given by ( |57| ) in a systematic manner when N is 
large and q > 1. We begin with the function B n (v,w;q) which is defined by the sum in 
equation (fSJ). The dominant terms are those around k = n/2, and so we may replace the 
g-binomial with the approximation (|69|) and also rewrite the sum as an integral over k. This 
gives for large n 

B n {v,w;q) ~ (-l)M(v, w; q) q* n *\vw\* n (CI) 

where A(v, w; q) is given by (ffl|). 

We now consider the sum (^) for Rn,u- In this summation we keep only the term with 
largest k as the others are exponentially suppressed. We find then that 



'N 



(q 1 vw; q x ) 



where we have expanded the product (W|V) using the identity 



{S N {v, w; q) + 2N3 N . 1 {v, w; q)) (C2) 



gS = (^')- (C3) 



which holds when q > 1 and for all x and where we have defined 

LfJ 



S N (v, w; q) = (l- q N ) £(-1) V^ W> «»-*-! B (l> „,. q) . (C4) 

{q;q)r{q;q)N-2r 



The main contribution to this latest sum SN{v,w;q) is where r is small. The approxi- 
mation 

(-^ vt't 1 *-gy^ <C5) 

{q;q)r{q;q)N-2r {q;q)r 

which follows from fl66|) and (|68|) is valid in that region and when combined with the asymp- 
totic expression for B n (v,w; q) yields 

LfJ 

S N {v, w; q) ~ (-l)^(t;, w; q) (1 - rtlH^V"" E T~T (~ ) " ( C6 ) 



r=0 



(g;g) r \mu/ 



We are now left with a single summation which may be estimated from the identity flC3|) . 
We see 



£t^(£J = (i/-; 9 -% + o( 9 -^) < ct > 



o (g;g) r \vw 
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and so to leading order in q we find 

S N (v, w; q) ~ (-l) N A(v, w; q) (l/vw; g" 1 ) «,M* A 'q^ (C8) 

Noting that SV-i is exponentially smaller than 5jv we may finally write down the asymp- 
totic form of Z N when q > 1 

^~A(^ ;g )(<rW/™;r%(^J ^ (C9) 



which is the expression presented in section 5.2 
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